{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.1.1 Preconditioners in NGSolve\n",
    "\n",
    "Preconditioners are approximate inverses which are used within iterative methods to solve linear or non-linear equations.\n",
    "\n",
    "Here are some built-in preconditioners in NGSolve:\n",
    "\n",
    "* Jacobi and Block-Jacobi\n",
    "* Direct solvers, i.e. sparse factorization\n",
    "* Multigrid with different block-smoothers\n",
    "* p-version element-level BDDC\n",
    "\n",
    "This tutorial quickly introduces how to use these within a solver. (In later tutorials, we will see how to measure condition numbers and customize preconditioners.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.csg import unit_cube\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "import netgen.gui\n",
    "%gui tk\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simple test problem \n",
    "\n",
    "We define a simple Poisson problem solver with the name of a preconditioner as argument in order to experiment with preconditioners.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveProblem(h=0.5, p=1, levels=1, \n",
    "                 condense = False,\n",
    "                 precond=\"local\"):\n",
    "    \"\"\"\n",
    "    Solve Poisson problem on l refinement levels.\n",
    "        h: coarse mesh size\n",
    "        p: polynomial degree \n",
    "        l: number of refinement levels\n",
    "        precond: name of a built-in preconditioner\n",
    "        condense: if true, perform static condensations\n",
    "    OUTPUT: \n",
    "        List of tuples of ndofs and iterations\n",
    "    \"\"\"\n",
    "    mesh = Mesh(unit_square.GenerateMesh(maxh=h))\n",
    "    # mesh = Mesh(unit_cube.GenerateMesh(maxh=h))\n",
    "    fes = H1(mesh, order=p, dirichlet=\"bottom|left\")\n",
    "    \n",
    "    u, v = fes.TnT() \n",
    "    a = BilinearForm(fes, eliminate_internal=condense)\n",
    "    a += SymbolicBFI(grad(u)*(grad(v)))\n",
    "    f = LinearForm(fes)\n",
    "    f += SymbolicLFI(1*v)\n",
    "    gfu = GridFunction(fes)\n",
    "    Draw (gfu)\n",
    "    c = Preconditioner(a, precond) # 'Register' c to a BEFORE assembly\n",
    "\n",
    "    steps = []\n",
    "    \n",
    "    for l in range(levels):\n",
    "        if l > 0: mesh.Refine()\n",
    "        fes.Update()\n",
    "        a.Assemble()\n",
    "        f.Assemble()\n",
    "        gfu.Update()\n",
    "\n",
    "        # Conjugate gradient solver\n",
    "        inv = CGSolver(a.mat, c.mat, maxsteps=1000)\n",
    "\n",
    "        # Solve steps depend on condense \n",
    "        if condense:\n",
    "            f.vec.data += a.harmonic_extension_trans * f.vec\n",
    "        gfu.vec.data = inv * f.vec\n",
    "        if condense:\n",
    "            gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "            gfu.vec.data += a.inner_solve * f.vec\n",
    "        steps.append ( (fes.ndof, inv.GetSteps()) )\n",
    "        Redraw ()\n",
    "    return steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Preconditioner` registers itself to the `BilinearForm`. Whenever the `BilinearForm` is re-assembled, the `Preconditioner` is updated as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `local` preconditioner \n",
    "\n",
    "The `local` preconditioner is a simple Jacobi preconditioner. \n",
    "The number of CG-iterations with the local preconditioner is proportional to $h^{-1} \\sim 2^l$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(8, 4)]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SolveProblem(precond=\"local\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(8, 4),\n",
       " (21, 9),\n",
       " (65, 25),\n",
       " (225, 51),\n",
       " (833, 103),\n",
       " (3201, 208),\n",
       " (12545, 418),\n",
       " (49665, 842)]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_local = SolveProblem(levels=8, precond=\"local\")\n",
    "res_local"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multigrid preconditioner \n",
    "\n",
    "A geometric multigrid `Preconditioner` uses the sequence of refined meshes to define a preconditioner of optimal iteration number (and complexity as well). It uses a direct solve on the coarsest level, and block Gauss-Seidel smoothers on the refined levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(8, 2),\n",
       " (21, 4),\n",
       " (65, 7),\n",
       " (225, 8),\n",
       " (833, 8),\n",
       " (3201, 8),\n",
       " (12545, 8),\n",
       " (49665, 8)]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res_mg = SolveProblem(levels=8, precond=\"multigrid\")\n",
    "res_mg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW9//H3ykwGwhTGEOaZiEAYFGVQURARnAW01avghLXtbYu29KpXq2h/rYrggIJcsYhaRbGKA/MoJAgoc0IIEKZMhCRAxrN+fxxERDAJOck+5+Tzeh6fsPfZ5+xvluHDztprr2WstYiIiP8KcLoAERGpXgp6ERE/p6AXEfFzCnoRET+noBcR8XMKehERP6egFxHxcwp6ERE/p6AXEfFzCnoRET8X5OkPNMZ0AR4BGgGLrbWvlveeRo0a2datW3u6FBERv7Zhw4Ysa21MeceZisx1Y4yZBVwHZFhru5+xfxjwEhAIvGmtnXLGawHA29baO8r7/ISEBJuUlFRuHSIi8iNjzAZrbUJ5x1W062Y2MOysEwQC04HhQFdgjDGm66nXrgc+Az6vRM0iIlINKhT01toVQM5Zu/sCKdbaVGttMTAPGHXq+AXW2uHAOE8WKyIilVeVPvoWwP4zttOBfsaYwcCNQCi/cEVvjJkATACIi4urQhkiIvJLPH4z1lq7DFhWgeNmADPA3Ud/9uslJSWkp6dTWFjo6RK9SlhYGLGxsQQHBztdioj4qaoE/QGg5Rnbsaf2eUR6ejpRUVG0bt0aY4ynPtarWGvJzs4mPT2dNm3aOF2OiNSwjLxCJr67kWlje9I4KqzazlOVcfSJQAdjTBtjTAhwO7CgMh9gjBlpjJlx7Nixn71WWFhIw4YN/TbkAYwxNGzY0O9/axGRc5u6OJnEtBymLkqu1vNU6IreGPMuMBhoZIxJBx631s40xkwEvsQ9vHKWtXZrZU5urf0U+DQhIWH8ec5bmY/zSbXhexSRn+o0eSFFpa7T2++s28c76/YRGhTAzqeHe/x8FQp6a+2Y8+z/HD8eQhkZGUlBQYHHPu+JJ54gMjKSP/zhDx77TBHxLYUlZfzmyvZMW7KbkyVlAIQFB3BNt6b8ZUSXajmnX02BkJFXyK2vryUjX10hIuJd8gtLeHXZbi57bgl//3IXEaGBGCA0KICiUhdRoUHV1k/vaND/Uh/9haiu/i5rLX/84x/p3r078fHxvPfee6dfe+6554iPj6dHjx48+uijALzxxhv06dOHHj16cNNNN3HixAmP1iMivuPo8WL++dVOBkxZwnNf7KBLs7rMm9Cf3q3qM65/K+Y/OIBx/VqRWVBUbTV4fHhlZZTXR/+DJz/dyraDeed9fX1aDmfO5PBDf5cx0Ld1g3O+p2vzujw+sluF6vzoo4/YtGkTmzdvJisriz59+jBw4EA2bdrEJ598wrp16wgPDycnx/1M2Y033sj48e5vafLkycycOZOHH364QucSEf9wJK+QN1akMnf9Pk4UlzGsW1MeHNKOi2LrAdC/bcPTxz49uvv5PsYjHA16T7k4th77ck5w9EQxLgsBBuqHhxDXINwjn79q1SrGjBlDYGAgTZo0YdCgQSQmJrJ8+XLuvvtuwsPd52nQwP2PypYtW5g8eTK5ubkUFBRwzTXXeKQOEfF++7JP8Ory3Xy4IZ0yaxnVozkPDG5HhyZRjtXkE0FfkSvvv8z/nrnr3Xeti8tcDO/elKdviK+B6n7urrvu4uOPP6ZHjx7Mnj2bZcuWOVKHiNScnYfzeXVZCgs2HyQoIIBbEmK5f1A7WnrogrMq/KaPPqugiHH9qqe/6/LLL+e9996jrKyMzMxMVqxYQd++fRk6dChvvfXW6T74H7pu8vPzadasGSUlJfzrX//yWB0i4n02789lwttJXPPiCr7adoR7L2/LqklD+NsN8V4R8uAjffQV8fqdP87U6en+rhtuuIG1a9fSo0cPjDE8//zzNG3alGHDhrFp0yYSEhIICQnh2muv5ZlnnuGpp56iX79+xMTE0K9fP/Lz8z1aj4g4y1rL2tRsXlm6m1UpWUTXCeaRKztw16WtqR8R4nR5P1Oh+eir27nmo9++fTtdulTPmFJvU5u+VxFfZq1lyY4Mpi9N4dt9uTSKDGX85W0Y178VkaE1f91c0fnofaKPXkTESWUuy2ffH+KVpSnsOJxPi3p1eGp0d27pHUtYcKDT5ZVLQS8ich7FpS7mb0zn1WW7Scs+QbuYCP5xSw+uv7g5wYG+87ypgl5E5Cwni8uYl7iPGStSOXSskO4t6vLaHb24umtTAgJ8b34qR4PeGDMSGNm+fXsnyxARASCvsIQ5a/cya9Ueso8X07dNA6bcdBEDOzTy6QkI/WbUjYjIhcouKGLW6j28vWYv+UWlDO4Uw0ND2tPnPE/W+xp13YhIrXUw9yRvrEzl3fX7KCp1cW33ZjwwuB3dW0Q7XZpHKehFpNbZk3Wc15bt5qON6VgLo3u24IHB7WgXE+l0adXC/4J+6bMw5DGnqxARL7T9UB6vLNvNZ98dJDgwgLF94xg/sC2x9b3jCdbq4jdTIJy2fIrHPiotLY3OnTtz11130bFjR8aNG8eiRYsYMGAAHTp0YP369WRmZjJ06FC6devGvffeS6tWrcjKyvJYDSJSdRv2HuWe2YkMf2klS3dkMGFgO1ZNuoInR3X3+5AHX7kZu/BROPx9xT/4rRHlH9M0HoaX/49CSkoKH3zwAbNmzaJPnz7MnTuXVatWsWDBAp555hliY2O54ooreOyxx/jiiy+YOXNmxesUkWpjrWV1SjbTl6awNjWbeuHB/H5oR359SWuiw4OdLq9G+UfXTe5eOLb/x+29q9xfo1tCvVZV+ug2bdoQH++eBbNbt25ceeWVGGOIj48nLS2NtLQ05s+fD8CwYcOoX79+lc4nIlXjclkWbT/C9GW72bw/lyZ1Q5k8ogtj+sYR4cA0Bd7AN77rClx5n/ZENDzhua6g0NDQ038OCAg4vR0QEEBpaSlBQb7RhCL+KCOvkInvbmTa2J40CA/hP98d4pVlKew6UkBcg3CeuSGem3q3IDTI+6cpqE5KqSoaMGAA77//PpMmTeKrr77i6NGjTpckUmv8sHzow3M3cuhYIftyTtCxSSQv3X4xI+KbEeRD0xRUJ/8L+kGP1ujpHn/8ccaMGcOcOXO45JJLaNq0KVFRzq0kI1IbdJq8kKJS1+ntdXvca0EEBxq+eGSgT05TUJ38L+g9OLSydevWbNmy5fT27Nmzf/ZaUVERX375JUFBQaxdu5bExMSfdPeIiGcdO1nCXZe2ZtbqPZSUuadZDwkMYHh8U/4yootC/hz8L+hr2L59+7j11ltxuVyEhITwxhtvOF2SiF/KKihi5qo9zFm7l4KiUlrUC+NgbiEhp5YPjQoNonFUmNNleiVNalZFHTp0YOPGjU6XIeK3DuaeZMaKVOYlnpqmIL4ZDw1uz0uLdzGkcxPG9o1j7vp9ZOYXOl2q1/KNcfQiUuukZR3n1V+YpqA6lw/1N17ddWOt9empQSvCG5ZyFPEmOw/nM31pCv/57iBBgQGM6RvHhFowTUF18tqgDwsLIzs7m4YNG/pt2Ftryc7OJixM/Yoim/bnMn1pCl9vO0JESCDjL2/LPZe3Ub+7B3ht0MfGxpKenk5mZqbTpVSrsLAwYmNjnS5DxBHWWr5JzWH60hRWpWQRXSeY317VgbsubU298BCny/MbXhv0wcHBtGnTxukyRKQaWGtZujOD6Ut3s2HvURpFhvLY8M6M69+KyFo6TUF1UouKSI0pc1m+2HKY6UtT2HYojxb16vC/o7pxa0JLwoJr9zQF1UlBLyLVrqTMxSebDvLKshRSM4/TtlEEf7/5Ikb3bEGwpimodgp6Eak2hSVlfJC0n9eWp3Ig9yRdmtVl2tieDO/ejEA9wVpj9MCUiHhcQVEpc9ft5Y2Ve8jML6JXXD2eGt2NIZ0a++0oOm9mvGEcd0JCgk1KSnK6DBGpotwTxcxek8Zbq9M4drKEAe0b8tCQ9lzS1n+HSTvJGLPBWptQ3nHquhGRKsvML+LNVam8s3Yvx4vLuKpLEx4a0o6ecVqIxxso6EXkgh3IPcnry3fzXuJ+SspcjLioOQ8ObkeXZnWdLk3OoKAXkUpLzSzg1WW7mb/xAAA39mrBA4Pb06ZRhMOVybko6EWkwrYdzGP6shQ+//4QIYEB3NG/FeMHtqVFvTpOlya/QEEvIuX6dt9Rpi9JYfGODCJDg7hvYDvuuawNMVFaZMcXKOhF5JystazZnc30pSms2Z1NvfBgfj+0I7++pDXR4cFOlyeVoKAXkZ+w1rJ4ewbTlqawaX8uMVGh/OXaLoztF0eE5qHxSfq/JlKLZeQVMvHdjUwb25OGEaF89v0hXlmawo7D+cTWr8PTo7tzc+9YzUPj4xT0IrXY1MXJJKbl8Ju5GzmSX8SerOO0i4ngH7f04PqLm2seGj+hoBephTpNXkhRqev09jd7cgAICjB89btBmofGz+ifa5Fa5mRxGfcPakto0I9//UMCAxh1cXPWPHaFQt4PaVIzkVrieFEpc77Zy5srU8kqKKZxVAiZ+cWEBAVQXOYiKjRIy/b5KUeD3lr7KfBpQkLCeCfrEPFneYUlvL0mjZmr9nD0RAmXd2jEb67swJsrU4mJCmNs3zjmrt9HZn6h06VKNVEfvYifOnaihFmr9/DW6j3kFZZyRefGPHxF+9MTjfVp3eD0sU+P7u5UmVIDFPQifibneDEzV6Xyf2v2UlBUytVdm/DwFR2Ij412ujRxiIJexE9k5hfx5spU5nyzl5MlZVzbvRkTr2ivmSRFQS/i647kFfL68lTmrt9LcamLkT2aM3FIezo0iXK6NPESCnoRH3Uw9ySvLd/NvMT9lLksoy9uwUND2tE2JtLp0sTLKOhFfMz+nBO8smw3/96wH2vh5t6xPDi4PXENw50uTbyUgl7ER6RlHWf60hQ+2niAQGO4rU9L7h/Ujtj6Cnj5ZQp6ES+XkpHPtCUpLNh8kODAAO7s34r7B7WjabQebpKKUdCLeKmdh/N5eUkyn31/iLCgQO65rA3jB7bV06tSaQp6ES+z5cAxXl6SzJdbjxAREsj9g9px72VtaBip1ZzkwijoRbzEpv25vLw4mcU7MogKC+I3V7Tnvy5rQ73wEKdLEx+noBdx2Ia9Oby0OIUVuzKJrhPMfw/tyK8ubU10HS3XJ56hoBdxyDep2UxdnMya3dk0iAhh0rDO3HlJKyK1XJ94mH6iRGqQtZbVKe6AX5+WQ6PIUCaPcK/HGh6iv45SPfSTJVIDrLUs25XJ1MXJbNyXS9O6YTwxsiu3943TeqxS7RT0ItXIWsui7Rm8vCSZ79KP0aKee8HtWxJiCQ1SwEvNUNCLVAOXy/LF1sO8vCSF7YfyiGsQznM3xXNDz1hCgrSCp9Qsjwe9MWY0MAKoC8y01n7l6XOIeJOMvEImvruRaWN70jAilP98d5DpS1PYdaSAto0i+MctPRh1cXOCAhXw4owKBb0xZhZwHZBhre1+xv5hwEtAIPCmtXaKtfZj4GNjTH3g/wEKevFrUxcnk5iWw2/nbeJwXiGpmcfp0DiSl26/mOsuaq7FtsVxFb2inw1MA97+YYcxJhCYDgwF0oFEY8wCa+22U4dMPvW6iF/qNHkhRaWu09trdmcDEBRg+PK3AwlQwIuXqNDvktbaFUDOWbv7AinW2lRrbTEwDxhl3J4DFlprv/VsuSLeobTMxZ+u6UR4yI83VEMCDaN6NGfNY1co5MWrVKWPvgWw/4ztdKAf8DBwFRBtjGlvrX3tXG82xkwAJgDExcVVoQyRmlNa5uKTTQd5eUkyadknqB8ezMniMkKCAigucxEVFqRJx8TrePxmrLV2KjC1AsfNAGYAJCQkWE/XIeJJZS7Lgs0HeHlxCqlZx+narC4z7uzNh9+mExMVxti+ccxdv4/M/EKnSxX5maoE/QGg5Rnbsaf2ifiNMpflP98d5KXFyaRmHqdz0yheu6M3V3dtQkCA4epuTU8f+/To7r/wSSLOqUrQJwIdjDFtcAf87cBYj1Ql4jCXy/Kf7w8xdXEyKRkFdGoSxSvjejGsW1P1v4vPqejwyneBwUAjY0w68Li1dqYxZiLwJe7hlbOstVsrc3JjzEhgZPv27StXtUg1cbksn285xEuLkknOKKBD40imje3Jtd2bKeDFZxlrne8eT0hIsElJSU6XIbXYD0+yvrQomZ1H8mkXE8EjV3VkRHwzjYMXr2WM2WCtTSjvOE2BILWatZYvtx7hxUW72HE4n7YxEXrQSfyOo0GvrhtxirWWr7cd4cVFyWw7lEebRhG8cFsPru/RQgEvfkddN1KrWGtZvD2DFxfvYsuBPFo3DOfhKzpoLhrxSeq6ETmDtZalOzN4cZF7uuC4BuH8/eaLuKFnCwW8+D0Fvfi1Hxb8eHFRMpv35xJbvw7P33QRN/RqQbACXmoJBb34JWstK5KzeOHrXWzan0uLenWYcmM8N/WOVcBLraObseJXrLWsSnEH/Lf73AH/zA3x3NxbC35I7aWbseIXrLWs2Z3NC1/vImnvUZpFh/HQkPZask/8mm7GSq2xdnc2Lyzaxfo9OTStG8ZTo7pxa5+WCniRUxT04rPWpboD/pvUHJrUDeXJ67txW5+WhAUr4EXOpKAXn5OYlsMLX+9ize5sYqJCeXxkV8b0jVPAi5yHbsaKz9iwN4cXvk5mVUoWjSJD+et1XRnXTwEvUh7djBWv9+2+o7zw9S5WJmfRMCKE+we1447+ragTooCX2k03Y8XnZOQVMvHdjUwb25PGUWFs2p/LC1/vYvmuTBpEhPDY8M7ceUkrwkP0YytSGfobI15j6uJkEtNy+J9PtlJUUsbSnZnUDw9m0rDO/OqSVkSE6sdV5ELob444rtPkhRSVuk5vf7HlMACBAYaVk64gUgEvUiV6VFAcN29Cf2Lr1zm9HRRgGBHflLWPKeRFPEGjbsQxOceLmb40hTlr91Lqcl/RhwQFUFLmon54CI2jwhyuUMQ/OBr01tpPgU8TEhLGO1mH1KwTxaW8tTqN15bt5nhxKbf0bsnhvJO0bBDB2L5xzF2/j8z8QqfLFPEb+r1YakxpmYv3k9J5cdEuMvKLGNq1CX+6phMdmkT95LinR3d3qEIR/6Sgl2rnXpf1MM9/sZPUrOMktKrPK+N6kdC6gdOlidQKCnqpVt+kZjNl4Q427c+lQ+NI3vhVAld1aYwxWpdVpKYo6KVabD+Ux/Nf7GDpzkya1g3j+Zsu4sZeWrZPxAkKevGo9KMn+OfXu5i/8QBRoUE8Orwzd13aWvPRiDhIwyvFI46eGir59tq9YGDCwLY8OKg90eHBTpcmUutpeKVUycniMmat3nN6qOTNvWP57VUdaV6vTvlvFpEaoa4buSClZS4+2OAeKnkkr4irujThT8M60fGsoZIi4jwFvVSKe6jkEZ7/cgepmcfp3ao+08b2oo+GSop4LQW9VNj6PTk8u3A7G/fl0i4mghl39mZo1yYaKini5RT0Uq6dh/N5/osdLN6RQZO6oTx3Uzw39YrVUEkRH6Ggl/M6kHuSF77exYffphMZGsSkYe6hklrZScS3KOjlZ3JPFPPKst3MXpMGwPjL2/Lg4HbUCw9xtjARuSAKejmtsKSMt1an8cqyFAqKSrmpVyy/G9qRFhoqKeLTFPRCaZmLD79N54WvkzmcV8iVnRvzx2Gd6Ny0rtOliYgH6MnYWsxay9fbjvD8lztJySigZ1w9Xrr9Yvq1beh0aSLiQcZa63QNJCQk2KSkJKfLqFUS03KYsnAHG/YepW1MBH+6pjPXdNNQSRFfYozZYK1NKO84dd34uYy8Qia+u5FpY3vSOCqMXUfyef6LnSzafoQmdUN59sZ4bumtoZIi/kxB7+emLk4mMS2HZz/fTnBgAP/ekE5EaBB/GtaJuy9to6GSIrWAgt5PdZq8kKJS1+nt+RsPAhBoDCv+OIT6ERoqKVJb6Pd1P7XyT0O4snNjfuhxDzQwtGtj1v75CoW8SC2joPdTqVnHWZmciQWCAw0uoElUGI2jwpwuTURqmLpu/NAHSfv58/zvCQkM4Nr4JkwY2I656/eRmV/odGki4gAFvR9xuSx//2onry7bzeUdGjFtbC+i67hXeHp6dHeHqxMRpyjo/cTJ4jJ+//4mFm45zNh+cTx5fTeCNWRSRFDQ+4WMvELGv53EdweOMXlEF+65rI0efBKR0xT0Pm7bwTzu/b9Eck+WMOPOBIZ2beJ0SSLiZRT0PmzJjiM8PHcjUWHBfHD/JXRrHu10SSLihRztxDXGjDTGzDh27JiTZfgcay2zVu3h3v9Lom1MJJ9MHKCQF5HzcjTorbWfWmsnREcrpCqqtMzF/3yylf/9zzaGdm3Ce/f1p0ldjY0XkfNT140PySssYeLcjazYlcl9g9oy6ZrOBATopquI/DIFvY/Yn3OC/5qdyJ6s4zx3Uzy39YlzuiQR8REKeh+wYe9R7puTREmZ5e17+nJpu0ZOlyQiPkRB7+UWbD7IHz7YTLPoMGbd1Yd2MZFOlyQiPkZB76WstUxdnMILi3bRt00DXr+jt2adFJELoqD3QoUlZTz64Xd8vOkgN/WK5ZkbuxMapAVCROTCKOi9THZBEffN2UDS3qP88ZpOPDi4naYzEJEqUdB7kZSMfO6enUhGXhHTx/ZixEXNnC5JRPyAgt5LrErO4oF/bSA0KJB5E/rTM66+0yWJiJ9Q0HuBuev28ddPttChcSRv/jqB2PrhTpckIn5EQe+gMpfl2c+38+aqPQzuFMPLY3oSFRbsdFki4mcU9A45XlTKI/M2sWj7Ee66tDWTR3QhSAuFiEg1UNA74NCxk9wzO4kdh/N48vpu/PrS1k6XJCJ+TEFfw75PP8a9bydyvKiMmXf1YUinxk6XJCJ+TkFfg77cepjfzttEg4gQPnygH52aRjldkojUAgr6GmCtZcaKVKZ8sYMesfV441cJxESFOl2WiNQSCvpqVlzq4q8fb+G9pP2MuKgZ/7ilB2HBms5ARGqOx4PeGNMW+AsQba292dOf70uOnSjh/nc2sDY1m4evaM/vruqohUJEpMZVaDyfMWaWMSbDGLPlrP3DjDE7jTEpxphHAay1qdbae6qjWF+SlnWcG15dzYa9R/nnrT3476s7KeRFxBEVHbg9Gxh25g5jTCAwHRgOdAXGGGO6erQ6H7UuNZvRr6zm6PFi3rm3Hzf2inW6JBGpxSoU9NbaFUDOWbv7AimnruCLgXnAqIqe2BgzwRiTZIxJyszMrHDB3u7DDencMXMdDSJC+PihAfRt08DpkkSklqvKo5gtgP1nbKcDLYwxDY0xrwE9jTGPne/N1toZ1toEa21CTExMFcpwXkZeIbe+toYnF2zlvz/YTJ/WDZj/wABaNYxwujQREc/fjLXWZgP3e/pzvdk/v97F+rSjrE87yu19WvLU6O4EazoDEfESVQn6A0DLM7ZjT+2rNTpNXkhRqesn++Yl7mf+xgPsfHq4Q1WJiPxUVS47E4EOxpg2xpgQ4HZgQWU+wBgz0hgz49ixY1Uowzkv3HYxIUE/jqQJCw5g1MXNWTlpiINViYj8VEWHV74LrAU6GWPSjTH3WGtLgYnAl8B24H1r7dbKnNxa+6m1dkJ0dHRl63aUtZY3V6by8LsbCQ0KxAChQQEUlbqICg2icVSY0yX6rqXPOl3BhfHVusG3awfVXwEVHXUzxlrbzFobbK2NtdbOPLX/c2ttR2ttO2vt36q3VO9woriU38zbxNOfbeeqLo3p27oB4/q3Yv6DAxjXrxWZBUVOl+jblk9xuoIL46t1g2/XDqq/Aoy1ttpPUp6EhASblJTkdBnlSss6zn1zNpCckc8frunEA4O0cHeFuVxwIhvyD0H+4Z9+LTjy023jg1NE2DLfrBt8u3bwj/qfuLDua2PMBmttQnnHOTrXjTFmJDCyffv2TpZRIUt2HOGReZsIDDDMvrsvAzv69pBQj7EWTuT8GNQFZ4b4GX8uOAKu0p+/P7whmAA4fsazFLbM/TXuEmg1oGa+jwuxdzXsW/vjtq/UDb5dO/hf/U+c6r4e9CgMOe+o9AumK/pyuFyWqUuSeXFRMt2a1+W1O3rTsoGPrem69NnK//BYCyeP/vxq+/TXIz8Ge1nxz99fpz5ENYOopu6vkU1+uh3V1L0vKOSn73si+oKvbhzlq3WDb9cOtbp+n7ii93bHTpbwu/c2sWRHBjf2asEzN8T75syTy6f8GPTWQlHeeYL7rEAvO8f9hrBoiGzqDupWl/40uH/4L7IpBOuGtIi3UNCfx47Dedw3ZwMHjp7kqVHduKN/K9/rjy8tgu8/cP951vAf+8NLTvz82JCoH4O6Zb+fBvfpK/CmEFLNv80MerR6P7+6+Grd4Nu1g+qvAEe7bs7oox+fnJzsWB1nW7D5IJP+/R1RYUG8ekcverfysflqCvPg/V9B6tKfvxbbB7qO+ml4RzWBUK12JeJrfKLrxlr7KfBpQkLCeCfr+EFJmYspC3cwc9Ue+rSuz/SxvWhc14e6IAoy4JtXIXEmFB2DNgNhwG/hnRt9uw9TRKpEXTenZOYXMXHut6zbk8Ndl7bmz9d2ISTIR+aryd4Na16GTXPdN0a7Xg8DHoEWvZ2uTES8gIIe2LjvKA+88y25J4t54bYe3NDTR+aPP7gRVr0I2xdAQBD0GAOX/gYanTVc1df7MEWkSmp10FtreXf9fp5YsJXGdUP58IFL6dbcy6djsNbd977qRdizHELrusO9/wPuPvdzqYZxuSLiO2rtA1OFJWU8/slW3kvaz8COMUy9/WLqhYeU/0anuMpg28ew+iU4tNl9E/WqJyHhbveQRxGR86iVD0wdyD3JA+9s4Lv0Y0wc0p7fDe1IoLeu51pyEjb9y90HfzQNGrZ3X8H3uB2CQp2uTkQc5BOjbpywJiWLie9upLjUxYw7e3N1t/N0dzjt5FFIfBPWve6eHqBFbxj6FHQeAQE++NCWiDim1gS9tZY3VqYyZeEO2sZE8vqdvWkXE+l0WT937AB88wpsmA3FBdD+KvcQydaXga89sCUiXqHASpSWAAAIM0lEQVRWBP3xolL+9O/v+Oz7Q1wb35Tnb+5BZKiXfeuZO2H1VPjuPbAu6H6je4hk03inKxMRH+dlaed5qZkF3DdnA7szC3hseGcmDGzrXVMZ7F/vHkGz8zMIquO+uXrJRKjfyunKRMRP+PWom6+3HeH3720iKNAw555+DGjfqFrOU2nWQvJX7oDft8Y90+OgSdB3AkR4SY0i4jf8ctRNmcvy4qJdvLwkhfgW0bx6Ry9i63vB1MJlJbDlQ/cQyYxtUDcWLp0IPe+EUC+8XyAiXq3WjrrJPVHMI/M2sXxXJrf0juWp0d2dn1q4+Dh8+zasnQ7H9kPjrnDD69D9JggMdrY2EfF7fhX02w7mcd87SRw+VsjfbujO2L5xzvbHH8+G9a/D+hnu4ZJxl8KIf0CHqzWCRkRqjE8HfUZeIRPf3ci0sT1ZnZLFYx99T3SdYN677xJ6xdWv2WLOXMXp6F5YOw2+nQOlJ6HTte4hknH9arYmERF8POinLk4mMS2HO99cx84jBfRt04DpY3sRE+XAE6PLp0CX69z971s+cq+DetFtMOA3ENOp5usRETnFJ4O+0+SFFJW6Tm/vPFIAwOb9uc6E/N5Ti/y+dhmERLonGOv/IES3qPlaRETO4uiE68aYkcaYGceOVW5RjJV/GsL1Fzcn8FQ3d3CgYdTFzVk5aUg1VPkLFj3pXtj3rWE/7isucIe9Ql5EvISjQW+t/dRaOyE6unKzLzauG0ZUaBAuICTQUOqyRIUG0TiqBleDSlvlHiqJcT/gBO5VnJ44pmmBRcSr+MgSSj+XVVDEuH6t+PihyxjXrxWZBUU1c+LiE7BwEswe4e6Hv3shXPO3mjm3iMgF8Mk+eoDX7/zxGYGnR3evmZPu+wY+fgByUqHvfXDV4xAS4X5NqziJiJfy2aCvUSWFsPRpWDMN6rWEX3/qXnj7TOquEREvpaAvT/oG+Ph+yNoFve+Gq5+C0CinqxIRqTAF/fmUFsHy52DVCxDVDO74CNpf6XRVIiKVpqA/l0ObYf4DkLEVet4B1zyjdVlFxGcp6M9UVgIr/wEr/g7hjWDs+9DxGqerEhGpEgX9D45shfn3w+Hv3FMXDJsC4Q2crkpEpMr8euGRCikrhdUvwrIpUKce3PYOdBnpXD0iIh7mk0/GekzmTpg5FJY8BZ1HwIPrFPIi4ndqZ9eNq8y9CMiSp90PPN38lnsxbhERP1T7gj57t/vp1v3roPN1cN0LENnY6apERKpN7Ql6l8u90tOiJyAoBG58A+Jv0UpPIuL3akfQ5+yBTybC3lXuZfxGToW6zZyuSkSkRvh30FsLSbPgq79CQCCMmg4Xj9NVvIjUKv4b9Ln7YcFESF0GbYfA9S+7JyQTEall/C/orYWNc+CLP4N1uW+29r5bV/EiUmv5V9DnHYRPH4Hkr6D15TBqGtRv7XRVIiKO8tkVpk5b+qz7Kn7zPHilP+xZCcOfh18tUMiLiOAPV/TLp8Dh72HnZ9CyP4x+BRq2c7oqERGv4dtBv3W++2vKIrj6aej/oHt0jYiInOabk5otfdZ9Jf+DsiL4ajIUFWhJPxGRsxhrrdM1kJCQYJOSkir/xrJSeKohPHHM80WJiHg5Y8wGa21Cecf59s3YQN/ueRIRqQm+HfQAgx51ugIREa/m+0GvPnkRkV/k+0EvIiK/SEEvIuLnFPQiIn5OQS8i4ucU9CIifs4rHpgyxmQCe8/aHQ2c60mos/eXt90IyPJAmedzvjo99b7yjqtoO51v/7mOO3Ofr7dfecdW9rWK7KvJ9jtfTZ58nyd/BivbfuD7P4PV2X6trLUx5VZgrfXK/4AZFdlfge0kJ+r01PvKO66i7VTR9jp7n6+3X3nHVva1iuyryfariTb05M9gZduvJtrQ19uvIv95c9fNpxXcX952dbvQ81X0feUdV9F2Ot/+cx1Xk21Y3e1X3rGVfa0i+/QzeP79ar+Kv+6xv5te0XVTnYwxSbYCc0HIuan9qkbtV3Vqw6rz5it6T5nhdAE+Tu1XNWq/qlMbVpHfX9GLiNR2teGKXkSkVlPQi4j4OQW9iIifq3VBb4xpa4yZaYz5t9O1+CJjzGhjzBvGmPeMMVc7XY+vMcZ0Mca8Zoz5tzHmAafr8UXGmAhjTJIx5jqna/EVfhH0xphZxpgMY8yWs/YPM8bsNMakGGMeBbDWplpr73GmUu9Uyfb72Fo7HrgfuM2Jer1NJdtvu7X2fuBWYIAT9XqbyrTfKZOA92u2St/mF0EPzAaGnbnDGBMITAeGA12BMcaYrjVfmk+YTeXbb/Kp16WS7WeMuR74DPi8Zsv0WrOpYPsZY4YC24CMmi7Sl/lF0FtrVwA5Z+3uC6ScuoIvBuYBo2q8OB9QmfYzbs8BC62139Z0rd6osj9/1toF1trhwLiardQ7VbL9BgP9gbHAeGOMX2RYdfPn1bVbAPvP2E4H+hljGgJ/A3oaYx6z1j7rSHXe75ztBzwMXAVEG2PaW2tfc6I4H3C+n7/BwI1AKLqi/yXnbD9r7UQAY8xdQJa11uVAbT7Hn4P+nKy12bj7l+UCWGunAlOdrsNXWWuXAcscLsPnWWtnO12DL/HnX3sOAC3P2I49tU8qRu1XNWq/qlH7eZA/B30i0MEY08YYEwLcDixwuCZfovarGrVf1aj9PMgvgt4Y8y6wFuhkjEk3xtxjrS0FJgJfAtuB9621W52s01up/apG7Vc1ar/qp0nNRET8nF9c0YuIyPkp6EVE/JyCXkTEzynoRUT8nIJeRMTPKehFRPycgl5ExM8p6EVE/JyCXkTEz/1/YiRLR2B+6yMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.plot(*zip(*res_local), \"-*\")\n",
    "plt.plot(*zip(*res_mg), \"-+\")\n",
    "plt.legend(['local', 'mg'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Multigrid implementation for higher order spaces \n",
    "\n",
    "For high order elements we use hierarchical spaces, where the (small) sub-spaces $V_E$, $V_F$, and $V_C$ are generated by basis functions associated with edges, faces, and cells:\n",
    "\n",
    "$$\n",
    "V_{hp} = V_{p=1} + \\sum_{\\text{edges }E} V_E + \\sum_{\\text{faces }F} V_F + \\sum_{\\text{cells }C} V_C\n",
    "$$\n",
    "\n",
    "The system matrix is stored as\n",
    "\n",
    "$$\n",
    "A = \\left( \\begin{array}{cccc}\n",
    "A_{VV} & A_{VE} & A_{VF} & A_{VC} \\\\\n",
    "A_{EV} & A_{EE} & A_{EF} & A_{EC} \\\\\n",
    "A_{FV} & A_{FE} & A_{FF} & A_{FC} \\\\\n",
    "A_{CV} & A_{CE} & A_{CF} & A_{CC} \\\\\n",
    "\\end{array} \\right)\n",
    "$$\n",
    "\n",
    "The $A_{VV}$-block is exactly the system matrix of a lowest order method. \n",
    "\n",
    "NGSolve's *multigrid implementation for a high order method uses h-version multigrid for the lowest order block,* and  local block-smoothing for the high order bubble basis functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 1 , res = [(8, 2), (21, 4), (65, 7), (225, 8)]\n",
      "p = 2 , res = [(21, 5), (65, 6), (225, 8), (833, 8)]\n",
      "p = 3 , res = [(40, 9), (133, 12), (481, 12), (1825, 13)]\n",
      "p = 4 , res = [(65, 12), (225, 15), (833, 16), (3201, 16)]\n",
      "p = 5 , res = [(96, 14), (341, 19), (1281, 20), (4961, 20)]\n",
      "p = 6 , res = [(133, 16), (481, 23), (1825, 23), (7105, 23)]\n",
      "p = 7 , res = [(176, 18), (645, 25), (2465, 26), (9633, 26)]\n",
      "p = 8 , res = [(225, 19), (833, 27), (3201, 28), (12545, 28)]\n",
      "p = 9 , res = [(280, 20), (1045, 29), (4033, 30), (15841, 30)]\n"
     ]
    }
   ],
   "source": [
    "for p in range(1,10):\n",
    "    r = SolveProblem(h=0.5, p=p, levels=4, condense=False, \n",
    "                     precond=\"multigrid\")\n",
    "    print (\"p =\",p,\", res =\",r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the number of iterations grows mildly with the order, and remains bounded with mesh refinement.\n",
    "\n",
    "Performing static condensation improves the situation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 1 , res = [(8, 2), (21, 4), (65, 7), (225, 8)]\n",
      "p = 2 , res = [(21, 5), (65, 6), (225, 8), (833, 8)]\n",
      "p = 3 , res = [(40, 5), (133, 6), (481, 7), (1825, 8)]\n",
      "p = 4 , res = [(65, 5), (225, 6), (833, 7), (3201, 8)]\n",
      "p = 5 , res = [(96, 5), (341, 6), (1281, 7), (4961, 8)]\n",
      "p = 6 , res = [(133, 5), (481, 6), (1825, 7), (7105, 8)]\n",
      "p = 7 , res = [(176, 5), (645, 6), (2465, 7), (9633, 8)]\n",
      "p = 8 , res = [(225, 5), (833, 6), (3201, 7), (12545, 8)]\n",
      "p = 9 , res = [(280, 5), (1045, 6), (4033, 7), (15841, 8)]\n"
     ]
    }
   ],
   "source": [
    "for p in range(1,10):\n",
    "    r = SolveProblem(h=0.5, p=p, levels=4, condense=True, \n",
    "                     precond=\"multigrid\")\n",
    "    print (\"p =\",p,\", res =\",r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Element-wise BDDC preconditioner\n",
    "\n",
    "A built-in element-wise BDDC (**B**alancing **D**omain **D**ecomposition preconditioner with **C**onstraints) preconditioner is also available. In contrast to local or multigrid preconditioners, the BDDC preconditioner needs access to the element matrices. This is exactly the reason why we register the preconditioner with the `BilinerForm`, and call the `bfa.Assemble()` later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 1 , res = [(27, 2), (89, 2), (321, 2)]\n",
      "p = 2 , res = [(89, 14), (321, 19), (1217, 20)]\n",
      "p = 3 , res = [(187, 18), (697, 23), (2689, 24)]\n",
      "p = 4 , res = [(321, 21), (1217, 28), (4737, 29)]\n",
      "p = 5 , res = [(491, 23), (1881, 31), (7361, 32)]\n",
      "p = 6 , res = [(697, 25), (2689, 34), (10561, 35)]\n",
      "p = 7 , res = [(939, 27), (3641, 36), (14337, 37)]\n",
      "p = 8 , res = [(1217, 28), (4737, 38), (18689, 40)]\n",
      "p = 9 , res = [(1531, 29), (5977, 39), (23617, 41)]\n"
     ]
    }
   ],
   "source": [
    "for p in range(1,10):\n",
    "    r = SolveProblem(h=0.25, p=p, levels=3, condense=True, \n",
    "                     precond=\"bddc\")\n",
    "    print (\"p =\",p,\", res =\",r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The BDDC preconditioner needs more iterations, but the work per iteration is less, so performance is similar to multigrid. **This element-wise BDDC preconditioner works well in shared memory parallel as well as in distributed memory mode.** Go to $\\S$[2.1.3](../unit-2.1.3-bddc/bddc.ipynb) for more about BDDC preconditioner. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
